## Separate Mechanisms (Results are estimated by lfe package version 2.8-7)
## These results are reported in Table 7
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)

##############################################################################################################################
########################
## Baseline Model
data_base1=data4
data_base2=data5

##############################################################################################################################
## Test seperate mechanisms for long sample
data_base1$CO2=data_base1$CO2_MASS..tons./1000
Cross_eff=felm(k_MWH~ factor(shutdown)  + ppt + tmean + tdmean + Wind  + SLOAD..1000.lbs.+ heat + CO2  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base1)
summary(Cross_eff)
Cross_CO2=felm(CO2~ factor(shutdown)   + ppt + tmean + tdmean + Wind + heat  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base1)
summary(Cross_CO2)
Cross_SO2=felm(SO2_MASS..tons.~ factor(shutdown)   + ppt + tmean + tdmean + Wind + heat + CO2 +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base1)
summary(Cross_SO2)
Cross_NOX=felm(NOX_MASS..tons.~ factor(shutdown)   + ppt + tmean + tdmean + Wind + heat + CO2   +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base1)
summary(Cross_NOX)
Cross_AOD=felm(AOD~ factor(shutdown)  + ppt + tmean + tdmean + Wind + heat + CO2  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base1)
summary(Cross_AOD)
stargazer(Cross_eff,Cross_CO2,Cross_SO2,Cross_NOX,Cross_AOD,digits=3)

##############################################################################################################################
## Test seperate mechanisms for short sample
data_base2$CO2=data_base2$CO2_MASS..tons./1000
Cross_eff=felm(k_MWH~ shutdown  + ppt + tmean + tdmean + Wind  + SLOAD..1000.lbs.+ heat + CO2  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base2)
summary(Cross_eff)
Cross_CO2=felm(CO2~ shutdown + ppt + tmean + tdmean + Wind + heat  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base2)
summary(Cross_CO2)
Cross_SO2=felm(SO2_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + heat + CO2 +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base2)
summary(Cross_SO2)
Cross_NOX=felm(NOX_MASS..tons.~ shutdown + ppt + tmean + tdmean + Wind + heat + CO2   +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base2)
summary(Cross_NOX)
Cross_AOD=felm(AOD~ shutdown + ppt + tmean + tdmean + Wind + heat + CO2  +factor(Year)+factor(weekday) | factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
               data = data_base2)
summary(Cross_AOD)
stargazer(Cross_eff,Cross_CO2,Cross_SO2,Cross_NOX,Cross_AOD,digits=3)
